---
title: "Mass survey"
author: "Endre Borbáth, Swen Hutter"
date: "20.07.2025"
output:
  html_document:
    toc: yes
---

```{r, echo = F, message = F, include=FALSE}

# R setup

rm(list = ls()) 

library(foreign)
library(dplyr)
library(tidyr)
library(ggplot2)
library(stringr)
library(haven)
library(forcats)
library(lubridate)
library(texreg)
library(DescTools)
library(broom)
library(ordinal)
library(here)
set.seed(1991)

cntry_ext <- c("DE", "IT", "PL")

```


```{r, include = F, echo = F, message = F}

# Data management

pol_controls <- c("contact_org", "contact_org2", "contact_org_orig") #"lrscale", "polinterest", "trust_state"

socdem_controls <- c("age", "age2", "edu_lev2", "gndr", "income", "region")

covid_specific <- c("home_child", "home_care", "infection", "neg_income_chng")

factors <- c("edu_lev2", "gndr", "region", "help_volunteer_past",
             "help_volunteer2", "help_volunteer",
             "support_unknown", "support_no_volunteer", "support_neighbor",
             "support_famfriends", "income", 
             "contact_org", "contact_org2", "contact_org_orig", 
             "org_membership", "home_child", "home_care", "infection")

dependent_vars <- c("help_volunteer", "help_volunteer2", "help_volunteer_past",
                    "support_unknown", "support_no_volunteer", "support_neighbor",
                    "support_famfriends", "vlntr_time_change")

load(here("study2_survey.Rdata"))

```

# Individual panel data

## Appendix A, figure 1

```{r, include = F, echo = F, message = F, warning=FALSE}

load(here("study2_overtime_survey.RData"))

dat1 <- dat %>% 
  mutate_at(vars(starts_with("help_")), ~ as.numeric(.)) %>% 
  pivot_longer(cols=starts_with("help_"), names_to = "help_form", values_to = "value") %>%
  mutate(type=case_when(grepl("_w1", help_form) ~ "wave1",
                        grepl("_w2", help_form) ~ "wave2",
                        TRUE ~ as.character(NA))) %>% 
  mutate(type2=ifelse(grepl("_past_w1", help_form), "wave0", NA)) %>% 
  mutate(type=case_when(type2=="wave0" ~ "wave0", 
                        TRUE ~ type)) %>% 
  select(-type2) %>% 
  mutate(help_form=case_when(grepl("_w1", help_form) ~ str_remove(help_form, "_w1"),
                        grepl("_w2", help_form) ~ str_remove(help_form, "_w2"),
                        TRUE ~ as.character(NA))) %>% 
  mutate(help_form=case_when(grepl("_past", help_form) ~ str_remove(help_form, "_past"),
                             TRUE ~ help_form)) %>% 
  pivot_wider(id_cols = c("n", "help_form"), names_from = type, values_from = "value")

dat2 <- dat %>%
  rename(weight_w1=wgt_2_w1,
         weight_w2=wgt_2_w1_adj) %>%
  select(n, weight_w1, weight_w2)

dat <- left_join(dat1, dat2) %>% 
  mutate(across(.cols = c(wave0, wave1, wave2),
           .fns = ~case_when(.x == 1 ~ 0,
                             .x >=2 ~ 1,
                             TRUE ~ as.numeric(NA)),
           .names = "{.col}_2")) %>% 
  group_by(help_form) %>% 
  mutate_at(vars("wave0", "wave1", "wave0_2", "wave1_2"), ~ weighted.mean(., w=weight_w1, na.rm = TRUE)) %>% 
  mutate_at(vars("wave2", "wave2_2"), ~ weighted.mean(., w=weight_w2, na.rm = TRUE)) %>% 
  select(-weight_w1, -weight_w2, -n) %>% 
  unique(.)
  
  

rm(dat1, dat2)


```

```{r, include = F, echo = F, message = F, warning=FALSE}


theme_set(theme_classic())

# Define functions. Source: https://github.com/jkeirstead/r-slopegraph
tufte_sort <- function(df, x="year", y="value", group="group", method="tufte", min.space=0.05) {
    ## First rename the columns for consistency
    ids <- match(c(x, y, group), names(df))
    df <- df[,ids]
    names(df) <- c("x", "y", "group")

    ## Expand grid to ensure every combination has a defined value
    tmp <- expand.grid(x=unique(df$x), group=unique(df$group))
    tmp <- merge(df, tmp, all.y=TRUE)
    df <- mutate(tmp, y=ifelse(is.na(y), 0, y))
  
    ## Cast into a matrix shape and arrange by first column
    require(reshape2)
    tmp <- dcast(df, group ~ x, value.var="y")
    ord <- order(tmp[,2])
    tmp <- tmp[ord,]
    
    min.space <- min.space*diff(range(tmp[,-1]))
    yshift <- numeric(nrow(tmp))
    ## Start at "bottom" row
    ## Repeat for rest of the rows until you hit the top
    for (i in 2:nrow(tmp)) {
        ## Shift subsequent row up by equal space so gap between
        ## two entries is >= minimum
        mat <- as.matrix(tmp[(i-1):i, -1])
        d.min <- min(diff(mat))
        yshift[i] <- ifelse(d.min < min.space, min.space - d.min, 0)
    }

    
    tmp <- cbind(tmp, yshift=cumsum(yshift))

    scale <- 1
    tmp <- melt(tmp, id=c("group", "yshift"), variable.name="x", value.name="y")
    ## Store these gaps in a separate variable so that they can be scaled ypos = a*yshift + y

    tmp <- transform(tmp, ypos=y + scale*yshift)
    return(tmp)
   
}

plot_slopegraph <- function(df) {
    ylabs <- subset(df, x==head(x,1))$group
    yvals <- subset(df, x==head(x,1))$ypos
    fontSize <- 3
    gg <- ggplot(df,aes(x=x,y=ypos)) +
        geom_line(aes(group=group),colour="grey80") +
        geom_point(colour="white",size=8) +
        geom_text(aes(label=y), size=fontSize, 
                  # family="American Typewriter"
                  ) +
        scale_y_continuous(name="", breaks=yvals, labels=ylabs)
    return(gg)
}    

```

```{r, include = F, echo = F, message = F, warning=FALSE}

dat1 <- dat %>% 
  select(help_form, wave1, wave0, wave2) %>% 
  pivot_longer(cols=c(wave1, wave0, wave2), values_to = "average", names_to = "period") %>% 
  mutate(period=case_when(period=="wave1" ~ 1,
                          period=="wave0" ~ 0,
                          period=="wave2" ~ 2)) %>% 
   filter(help_form!="help_other") %>% 
  mutate(help_form=case_when(help_form=="help_shopping" ~ "Shopping",
help_form=="help_childcare" ~ "Childcare",
help_form=="help_general" ~ "General",
help_form=="help_emo" ~ "Emotional",
help_form=="help_symb" ~ "Symbolic",
help_form=="help_fin" ~ "Financial",
help_form=="help_donation" ~ "Donation",
help_form=="help_volunteer" ~ "Volunteering",
help_form=="help_other" ~ "Other"))

```

```{r, include = T, echo = T, message = F, warning=FALSE, fig.width=9, fig.height=8}

## Prepare data    
df <- tufte_sort(dat1, 
                 x="period", 
                 y="average", 
                 group="help_form", 
                 method="tufte", 
                 min.space=0.05)

df <- transform(df, 
                x=factor(x, levels=c(0, 1, 2), 
                            labels=c("Pre-crisis\nMarch 2019-\nMarch 2020","Wave 1\nOctober 14. 2020-\nNovember 04. 2020","Wave 2\nMarch 02. 2021-\nMarch 11. 2021")), 
                y=round(y, digits = 2))

## Plot
p <- plot_slopegraph(df) + 
  # labs(title="Average level of social support in the\nperiod of the COVID-19 crisis",
  #                               caption="Average level of social support offered by respondents to someone living\noutside of their household. Scale: 1-5.") + 
                      theme(axis.title=element_blank(),
                            axis.ticks = element_blank(),
                            axis.text = element_text(colour="black")) 

p

# ggsave(filename="help_average_support.png",
#        plot=p,
#        path=fig_path,
#        scale=0.8,
#        width=9,
#        height=8)

```

## Figure 2, main text

```{r, include = F, echo = F, message = F, warning=FALSE}

dat1 <- dat %>% 
  select(help_form, wave1_2, wave0_2, wave2_2) %>% 
  pivot_longer(cols=c(wave1_2, wave0_2, wave2_2), values_to = "average", names_to = "period") %>% 
  mutate(period=case_when(period=="wave1_2" ~ 1,
                          period=="wave0_2" ~ 0,
                          period=="wave2_2" ~ 2)) %>% 
  filter(help_form!="help_other") %>% 
  mutate(help_form=case_when(help_form=="help_shopping" ~ "Shopping",
help_form=="help_childcare" ~ "Childcare",
help_form=="help_general" ~ "General",
help_form=="help_emo" ~ "Emotional",
help_form=="help_symb" ~ "Symbolic",
help_form=="help_fin" ~ "Financial",
help_form=="help_donation" ~ "Donation",
help_form=="help_volunteer" ~ "Volunteering",
help_form=="help_other" ~ "Other"))

```


```{r, include = T, echo = T, message = F, warning=FALSE, fig.width=9, fig.height=8}

## Prepare data    
df <- tufte_sort(dat1, 
                 x="period", 
                 y="average", 
                 group="help_form", 
                 method="tufte", 
                 min.space=0.05)

df <- transform(df, 
                x=factor(x, levels=c(0, 1, 2), 
                            labels=c("Pre-crisis\nMarch 2019-\nMarch 2020","Wave 1\nOctober 14. 2020-\nNovember 04. 2020","Wave 2\nMarch 02. 2021-\nMarch 11. 2021")), 
                y=round(y, digits = 2))

## Plot
p <- plot_slopegraph(df) + 
  # labs(title="Proportion of those who offered social support\nin the period of the COVID-19 crisis",
  #                               caption="Weighted proportion of those who replied supporting 'rarely', 'sometimes',\n'often', or 'very often' someone living outside of their household.") + 
                      theme(axis.title=element_blank(),
                            axis.ticks = element_blank(),
                            axis.text = element_text(colour="black")) 

p

# ggsave(filename="help_prop_support.png",
#        plot=p,
#        path=fig_path,
#        scale=0.8,
#        width=9,
#        height=8)

```

```{r, include = F, echo = F, message = F, warning=FALSE}

level_contact_country <- reg_dat |> 
  select(cntry, contact_org, contact_org2, wgt_2) |> 
  mutate(contact_org2=ifelse(contact_org==1, wgt_2, 0)) |> 
  mutate(seldom=ifelse(contact_org==2, wgt_2, 0)) |> 
  mutate(sometimes=ifelse(contact_org==3, wgt_2, 0)) |> 
  mutate(often=ifelse(contact_org==4, wgt_2, 0)) |> 
  select(-contact_org) |> 
  group_by(cntry) |> 
  mutate_at(vars(contact_org2, seldom, sometimes, often, wgt_2), ~sum(., na.rm = TRUE)) |> 
  ungroup() |> 
  distinct() |> 
  mutate_at(vars(contact_org2, seldom, sometimes, often), ~.*100/wgt_2) |> 
  mutate(contact_org2=100-contact_org2) |> 
  select(-wgt_2)
  
level_contact <- level_contact_country |> 
  summarise_all(mean)

```



```{r, include = F, echo = F, message = F, warning=FALSE}

prop_table <- reg_dat %>% 
  mutate(contact=as.numeric(ifelse(contact_org2==1, wgt_2, 0))) %>% 
  select(cntry, contact, wgt_2) %>% 
  group_by(cntry) %>% 
  mutate_at(vars(wgt_2, contact), ~sum(., na.rm=TRUE)) %>% 
  ungroup(.) %>% 
  unique(.) %>% 
  mutate(contact_share=contact*100/wgt_2)

```


# Volunteering/ supporting others - reference: the whole sample

## With control variables

```{r, include = F, echo = F, message = F, warning=FALSE}

# Regression models

M1 <- glm(as.formula(paste0("help_volunteer2 ~ contact_org + org_membership + cntry + ",
                            paste(covid_specific, collapse = "+"), "+",
                            paste(socdem_controls, collapse = "+"))),
    data = reg_dat,
    weights = wgt_2,
    family = "binomial")

M2 <- glm(as.formula(paste0("support_unknown ~ contact_org + org_membership + cntry + ",
                            paste(covid_specific, collapse = "+"), "+",
                            paste(socdem_controls, collapse = "+"))),
    data = reg_dat,
    weights = wgt_2,
    family = "binomial")

M3 <- glm(as.formula(paste0("support_no_volunteer ~ contact_org + org_membership + cntry + ",
                            paste(covid_specific, collapse = "+"), "+",
                            paste(socdem_controls, collapse = "+"))),
    data = reg_dat,
    weights = wgt_2,
    family = "binomial")

M4 <- glm(as.formula(paste0("support_neighbor ~ contact_org + org_membership + cntry + ",
                            paste(covid_specific, collapse = "+"), "+",
                            paste(socdem_controls, collapse = "+"))),
    data = reg_dat,
    weights = wgt_2,
    family = "binomial")

M5 <- glm(as.formula(paste0("support_famfriends ~ contact_org + org_membership + cntry + ",
                            paste(covid_specific, collapse = "+"), "+",
                            paste(socdem_controls, collapse = "+"))),
    data = reg_dat,
    weights = wgt_2,
    family = "binomial")

```

# Table 3, appendix A

```{r, results='asis', include = T, echo = F, message = F, warning=FALSE}

R2s <- c(PseudoR2(M1, "McFadden"), PseudoR2(M2, "McFadden"),
         PseudoR2(M4, "McFadden"), PseudoR2(M5, "McFadden"))

htmlreg(l=list(M1, M2, M4, M5), doctype = FALSE, center = TRUE, siunitx=TRUE, 
       use.packages=FALSE, booktabs = TRUE, float.pos = "h!",
       caption = "Individual-level differences in solidarity action",
       single.row=TRUE, no.margin = TRUE,
       custom.model.names = c(" ", "strangers",
                              "neighbours", "family/ friends"),
       custom.header = list("Volunteering" = 1, "Helping" = 2,
                            "Helping" = 3, "Helping" = 4),
       custom.gof.rows = list("McFadden Sq(R)" =  R2s),
       custom.coef.map = list(#"(Intercept)" = "Intercept",
                              "contact_org2"="Seldom (ref: never)", 
                              "contact_org3"="Sometimes",
                              "contact_org4"="Often",
                              "infection1" = "Infection (self/ in household)",
                              "home_child1" = "Children at home",
                              "home_care1" = "Care work at home",
                              "neg_income_chng" = "Negative income change",
                              "org_membership1" = "Member in civ. soc. org",
                              "age" = "Age",
                              "age2" = "Sq(Age)",
                              "edu_lev21" = "High education",
                              "gndr2" = "Female",
                              "income2" = "Coping",
                              "income3" = "Difficult",
                              "income4" = "Very difficult",
                              "region2" = "Mid. town",
                              "region3" = "Small town",
                              "region4" = "Countryside",
                              "cntryIT" = "Italy",
                              "cntryPL" = "Poland"),
       groups = list("Contacted by an org." = 1:3,
                     "Crisis specific" = 4:7,
                     "Socio-demographics" = 8:12,
                     "Income (ref: comfortable)" = 13:15,
                     "Region (ref: big city)" = 16:18,
                     "Country FE (ref: Germany)" = 19:20),
       include.deviance = FALSE,
       include.loglik = FALSE,
       # file=paste0(fig_path, "general_participation.tex"),
      custom.note	= "%stars")

```

Results are based on the first wave, weighted by the socio-demographic weight. These are logit models. In the case of volunteering, I coded never as 0 and everything else as 1. In the helping questions I coded everyone who did not help as zero.

\newpage

Predicted probability of being contacted by an organization

```{r, include = F, echo = F, message = F, warning=FALSE}

new_dat <- data.frame("age" = mean(reg_dat$age),  
                      "age2" = mean(reg_dat$age2),
                      "home_child" = as.factor(0),
                      "home_care" = as.factor(0),
                      "infection" = as.factor(0),
                      "neg_income_chng" = mean(reg_dat$neg_income_chng),
                      "edu_lev2" = as.factor(0), 
                      "gndr" = as.factor(1), 
                      "income" = as.factor(1), 
                      "region" = as.factor(1),
                      "org_membership" = as.factor(0),
                      "cntry" = as.factor("DE"),
                      "contact_org" = as.factor(seq(1, 4, 1)))

results <- NULL

for (i in paste0("M", seq(1, 5, 1))) {

# Calculate the predicted probabilities

# probs <- as.data.frame(predict(!!(i), newdata = new_dat, type = "response", se.fit = TRUE))
probs <- as.data.frame(do.call("predict", list(as.name(i), newdata = new_dat, type = "response", se.fit = TRUE)))
probs$contact_org <- seq.int(nrow(probs)) 
probs$models <- i
probs$residual.scale <- NULL

results <- rbind(results, probs)

}

```


### Figure 4, main text

```{r, include = T, echo = F, message = F, warning=FALSE, fig.width=5, fig.height=5, fig.align="center"}

results <- results %>%
  filter(models!="M3") %>% 
  mutate(type=case_when(models=="M1" ~ "Volunteering",
                        models=="M2" ~ "Helping strangers",
                        models=="M4" ~ "Helping neighbours",
                        models=="M5" ~ "Helping family or friends")) %>% 
  mutate(type=factor(type, levels=c("Volunteering", "Helping strangers", 
                                   "Helping neighbours", "Helping family or friends"))) %>% 
  mutate(pm=fit,
         pu=fit + se.fit * 1.96, # 95% confidence interval
         pl=fit - se.fit * 1.96) %>% # 95% confidence interval 
  mutate(contact_org=as.factor(contact_org))

ggplot(results, aes(x=contact_org, y=pm, group=1)) +
  geom_point() + 
  geom_line(linetype="dashed") +
  geom_errorbar(aes(ymin=pl, ymax=pu), width=.2,
                 position=position_dodge(0.05)) +
  facet_wrap(~ type, ncol = 2) +
  scale_x_discrete(labels=c("1"="Never", "2"="Seldom", "3"="Sometimes",
                            "4"="Often")) +
  theme_linedraw() +
  ylab("Predicted probability") + ylim(0, 1) +
  xlab("Being contacted by an organization")
  # theme(axis.title.x=element_blank())

# ggsave(filename="margins_whole_sample.png",
#        plot=last_plot(),
#        path=fig_path,
#        scale=1.1,
#        width=5,
#        height=5)

```

# Table 4, Appendix A

```{r, results='asis', include = T, echo = F, message = F, warning=FALSE}

R2s <- c(PseudoR2(M1, "McFadden"), PseudoR2(M2, "McFadden"),
         PseudoR2(M3, "McFadden"))

htmlreg(l=list(M1, M2, M3), doctype = FALSE, center = TRUE, siunitx=TRUE, 
       use.packages=FALSE, booktabs = TRUE, float.pos = "h!",
       caption = "Individual-level differences in solidarity action",
       single.row=TRUE, no.margin = TRUE,
       custom.model.names = c(" ", " ",
                              "(without volunteering)"),
       custom.header = list("Volunteering" = 1, "Helping strangers" = 2,
                            "Helping strangers" = 3),
       custom.gof.rows = list("McFadden Sq(R)" =  R2s),
       custom.coef.map = list(#"(Intercept)" = "Intercept",
                              "contact_org2"="Seldom (ref: never)", 
                              "contact_org3"="Sometimes",
                              "contact_org4"="Often",
                              "infection1" = "Infection (self/ in household)",
                              "home_child1" = "Children at home",
                              "home_care1" = "Care work at home",
                              "neg_income_chng" = "Negative income change",
                              "org_membership1" = "Member in civ. soc. org",
                              "age" = "Age",
                              "age2" = "Sq(Age)",
                              "edu_lev21" = "High education",
                              "gndr2" = "Female",
                              "income2" = "Coping",
                              "income3" = "Difficult",
                              "income4" = "Very difficult",
                              "region2" = "Mid. town",
                              "region3" = "Small town",
                              "region4" = "Countryside",
                              "cntryIT" = "Italy",
                              "cntryPL" = "Poland"),
       groups = list("Contacted by an org." = 1:3,
                     "Crisis specific" = 4:7,
                     "Socio-demographics" = 8:12,
                     "Income (ref: comfortable)" = 13:15,
                     "Region (ref: big city)" = 16:18,
                     "Country FE (ref: Germany)" = 19:20),
       include.deviance = FALSE,
       include.loglik = FALSE,
       # file=paste0(fig_path, "general_participation.tex"),
      custom.note	= "%stars")

```

Results are based on the first wave, weighted by the socio-demographic weight. These are logit models. In the case of volunteering, I coded never as 0 and everything else as 1. In the helping questions I coded everyone who did not help as zero.

\newpage

Predicted probability of being contacted by an organization

```{r, include = F, echo = F, message = F, warning=FALSE}

new_dat <- data.frame("age" = mean(reg_dat$age),  
                      "age2" = mean(reg_dat$age2),
                      "home_child" = as.factor(0),
                      "home_care" = as.factor(0),
                      "infection" = as.factor(0),
                      "neg_income_chng" = mean(reg_dat$neg_income_chng),
                      "edu_lev2" = as.factor(0), 
                      "gndr" = as.factor(1), 
                      "income" = as.factor(1), 
                      "region" = as.factor(1),
                      "org_membership" = as.factor(0),
                      "cntry" = as.factor("DE"),
                      "contact_org" = as.factor(seq(1, 4, 1)))

results <- NULL

for (i in paste0("M", seq(1, 3, 1))) {

# Calculate the predicted probabilities

# probs <- as.data.frame(predict(!!(i), newdata = new_dat, type = "response", se.fit = TRUE))
probs <- as.data.frame(do.call("predict", list(as.name(i), newdata = new_dat, type = "response", se.fit = TRUE)))
probs$contact_org <- seq.int(nrow(probs)) 
probs$models <- i
probs$residual.scale <- NULL

results <- rbind(results, probs)

}

```

### Figure 3, appendix A

```{r, include = T, echo = F, message = F, warning=FALSE, fig.width=5, fig.height=5, fig.align="center"}

results <- results %>%
  mutate(type=case_when(models=="M1" ~ "Volunteering",
                        models=="M2" ~ "Helping strangers",
                        models=="M3" ~ "Helping strangers (without volunteering)")) %>% 
  mutate(type=factor(type, levels=c("Volunteering", "Helping strangers", 
                                    "Helping strangers (without volunteering)"))) %>% 
  mutate(pm=fit,
         pu=fit + se.fit * 1.96, # 95% confidence interval
         pl=fit - se.fit * 1.96) %>% # 95% confidence interval 
  mutate(contact_org=as.factor(contact_org))

ggplot(results, aes(x=contact_org, y=pm, group=1)) +
  geom_point() + 
  geom_line(linetype="dashed") +
  geom_errorbar(aes(ymin=pl, ymax=pu), width=.2,
                 position=position_dodge(0.05)) +
  facet_wrap(~ type, ncol = 3) +
  scale_x_discrete(labels=c("1"="Never", "2"="Seldom", "3"="Sometimes",
                            "4"="Often")) +
  theme_linedraw() +
  ylab("Predicted probability") + ylim(0, 1) +
  xlab("Being contacted by an organization")
  # theme(axis.title.x=element_blank())

# ggsave(filename="margins_whole_sample_appendix.png",
#        plot=last_plot(),
#        path=fig_path,
#        width=6.5,
#        height=3,
#        scale=1.2,
#        dpi=400)

```


# Being contacted and change in the level of volunteering as DVs

## Being contacted by an organization as a DV

```{r, include = F, echo = F, message = F, warning=FALSE}

M1 <- lm(as.formula(paste0("as.numeric(contact_org_orig) ~ cntry + ",
                            paste(covid_specific, collapse = "+"), "+",
                            paste(socdem_controls, collapse = "+"))),
    data = subset(reg_dat, org_membership==0),
    weights = wgt_2)

M2 <- lm(as.formula(paste0("as.numeric(contact_org_orig) ~ cntry + ",
                            paste(covid_specific, collapse = "+"), "+",
                            paste(socdem_controls, collapse = "+"))),
    data = subset(reg_dat, org_membership==1),
    weights = wgt_2)

```

## Table 1, Appendix A

```{r, results='asis', include = T, echo = F, message = F, warning=FALSE}

htmlreg(l=list(M1, M2), doctype = FALSE, center = TRUE, siunitx=TRUE, 
       use.packages=FALSE, booktabs = TRUE, float.pos = "h!",
       caption = "Individual-level differences in being contacted",
       single.row=TRUE, no.margin = TRUE,
       custom.model.names = c("Non members", "Members"),
       custom.coef.map = list("infection1" = "Infection (self/ in household)",
                              "home_child1" = "Children at home",
                              "home_care1" = "Care work at home",
                              "neg_income_chng" = "Negative income change",
                              "age" = "Age",
                              "age2" = "Sq(Age)",
                              "edu_lev21" = "High education",
                              "gndr2" = "Female",
                              "income2" = "Coping",
                              "income3" = "Difficult",
                              "income4" = "Very difficult",
                              "region2" = "Mid. town",
                              "region3" = "Small town",
                              "region4" = "Countryside",
                              "cntryIT" = "Italy",
                              "cntryPL" = "Poland"),
       groups = list("Crisis specific" = 1:4,
                     "Socio-demographics" = 5:8,
                     "Income (ref: comfortable)" = 9:11,
                     "Region (ref: big city)" = 12:14,
                     "Country FE (ref: Germany)" = 15:16),
       include.deviance = FALSE,
       include.loglik = FALSE,
       # file=paste0(fig_path, "general_participation.tex"),
      custom.note	= "%stars")

```

Results are based on the first wave, weighted by the socio-demographic weight. 
An effect of 1 means a 1 point step up on the five points scale above, it's a different interpretation than in the case of probabilities in the logit model above. 
I include below the effect plot.

\newpage

## Figure 3, main text

```{r, include = F, echo = F, message = F, warning=FALSE}

mod1_res <- tidy(M1)
mod1_res <- mod1_res %>% 
  mutate(type="No member")
mod1_ci <- as.data.frame(confint(M1))
mod1_ci$term <- rownames(mod1_ci)
mod1_res <- merge(mod1_ci, mod1_res, by="term", all=TRUE)
mod1_ci <- as.data.frame(confint(M1, level=0.9))
mod1_ci$term <- rownames(mod1_ci)
mod1_res <- merge(mod1_ci, mod1_res, by="term", all=TRUE)

mod2_res <- tidy(M2)
mod2_res <- mod2_res %>% 
  mutate(type="Member")
mod2_ci <- as.data.frame(confint(M2))
mod2_ci$term <- rownames(mod2_ci)
mod2_res <- merge(mod2_ci, mod2_res, by="term", all=TRUE)
mod2_ci <- as.data.frame(confint(M2, level=0.9))
mod2_ci$term <- rownames(mod2_ci)
mod2_res <- merge(mod2_ci, mod2_res, by="term", all=TRUE)

probs <- bind_rows(mod1_res, mod2_res)

```


```{r, include = T, echo = F, message = F, warning=FALSE, fig.width=5, fig.height=5}

plot_dat <- probs %>% 
  filter(term!="(Intercept)") %>% 
  filter(!grepl("cntry", term)) %>% 
  mutate(term=case_when(term=="edu_lev21" ~ "High educ.",
                        term=="gndr2" ~ "Female",
                        term=="region2" ~ "Mid. town",
                        term=="region3" ~ "Small town",
                        term=="region4" ~ "Countryside",
                        term=="age" ~ "Age",
                        term=="age2" ~ "Sq(Age)",
                        term=="infection1" ~ "Infection (self/ in household)",
                        term=="home_child1" ~ "Children at home",
                        term=="home_care1" ~ "Care work at home",
                        term=="neg_income_chng" ~ "Negative income change",
                        term=="income2" ~ "Coping",
                        term=="income3" ~ "Difficult",
                        term=="income4" ~ "Very difficult",
                        TRUE ~ term)) 

titles <- plot_dat %>%
  filter(term %in% c("Age", "Coping", "Difficult", "Female")) %>%  # I just need some unique values
  mutate_at(vars(c(2:9)), ~ NA)

titles$term <- rep(c("Crisis specific", "Socio-demographics", "Income (ref: comfortable)", "Region (ref: big city)"), 1)

plot_dat <- bind_rows(plot_dat, titles) %>% 
  mutate(term=factor(term, levels = c("Crisis specific",
                                      "Infection (self/ in household)",
                                      "Children at home",
                                      "Care work at home",
                                      "Negative income change",
                                      "Socio-demographics",
                                      "Age", "Sq(Age)", "High educ.", 
                                      "Female", 
                                      "Income (ref: comfortable)",
                                      "Coping", "Difficult", "Very difficult",
                                      "Region (ref: big city)",
                                      "Mid. town", "Small town", "Countryside"))) %>%
  mutate(term=fct_rev(term)) %>% 
  mutate(new_labels = as.character(term)) 

space_between_bars <- 0.5

res1 <- ggplot(plot_dat, aes(y=term, color=type)) +
  geom_vline(xintercept = 0, linetype="dashed", color="gray60") +
  geom_point(aes(x=estimate, shape=type), position=position_dodge(width=space_between_bars), size=2) +
  geom_linerange(aes( xmin=`2.5 %`, xmax=`97.5 %`, 
                     group=type), size=0.5, position=position_dodge(width=space_between_bars), key_glyph = "path") +
  geom_linerange(aes(xmin=`5 %`, xmax=`95 %`,
                     group=type), size=1.1, position=position_dodge(width=space_between_bars), key_glyph = "path") +
  xlab("OLS Estimates") +
  scale_color_manual(name="", values=c("#758fc7", "#cc592d"), guide=guide_legend(reverse = TRUE, nrow=1, byrow=TRUE)) +
  scale_shape_discrete(name="", guide=guide_legend(reverse = TRUE, nrow=1, byrow=TRUE)) +
  scale_y_discrete(labels = function(y) stringr::str_wrap(y, width = 22)) +
  theme_linedraw() +
  theme(legend.title=element_blank(), 
        legend.position="bottom",
        legend.key.width = unit(1,"cm"),
        axis.title.y = element_blank(),
        axis.text.y=element_text(face=c(rep('plain', 3), 'bold', 
                                        rep('plain', 3), 'bold', 
                                        rep('plain', 4), 'bold', 
                                        rep('plain', 4), 'bold'),
                                 hjust=1),
        # axis.text.y = element_text(face = c(rep(2, 4))),
        axis.ticks.y = element_blank())

# ggsave(plot=last_plot(),
#        filename="contact_by_membership.png",
#        path = fig_path,
#        width=5,
#        height=5,
#        scale=3,
#        dpi=400,
#        units="cm")



```

\newpage

I also do now the regression to see how do volunteers and those who supported previously unknown people differ from the general population.. In this part I only work with volunteering and supporting unknown people as the dependent variable - again dichotomized. To get to a sample of the ones that were reached by organizations, I dichotomize the contact variable and distinguish those who report some form of contact (one third of the total sample). I then split the sample and run the same model both on the full sample and on the subset of those who were contacted.

### Table 2, Appendix A

\newpage

## Volunteering

```{r, include = F, echo = F, message = F, warning=FALSE}

# Regression models

M1 <- glm(as.formula(paste0("help_volunteer2 ~ org_membership + cntry + ",
                            paste(covid_specific, collapse = "+"), "+",
                            paste(socdem_controls, collapse = "+"))),
    data = subset(reg_dat, as.numeric(contact_org)==1),
    weights = wgt_2,
    family = "binomial")

M2 <- glm(as.formula(paste0("help_volunteer2 ~ org_membership + cntry + ",
                            paste(covid_specific, collapse = "+"), "+",
                            paste(socdem_controls, collapse = "+"))),
    data = subset(reg_dat, as.numeric(contact_org)>1),
    weights = wgt_2,
    family = "binomial")

```


```{r, results='asis', include = T, echo = F, message = F, warning=FALSE}

R2s <- c(PseudoR2(M1, "McFadden"), PseudoR2(M2, "McFadden"))

htmlreg(l=list(M1, M2), doctype = FALSE, center = TRUE, siunitx=TRUE, 
       use.packages=FALSE, booktabs = TRUE, float.pos = "h!",
       caption = "Individual-level differences in volunteering",
       single.row=TRUE, no.margin = TRUE,
       custom.model.names = c("Not contacted", "Contacted"),
       custom.gof.rows = list("McFadden Sq(R)" =  R2s),
       custom.coef.map = list("infection1" = "Infection (self/ in household)",
                              "home_child1" = "Children at home",
                              "home_care1" = "Care work at home",
                              "neg_income_chng" = "Negative income change",
                              "org_membership1" = "Member in civ. soc. org",
                              "age" = "Age",
                              "age2" = "Sq(Age)",
                              "edu_lev21" = "High education",
                              "gndr2" = "Female",
                              "income2" = "Coping",
                              "income3" = "Difficult",
                              "income4" = "Very difficult",
                              "region2" = "Mid. town",
                              "region3" = "Small town",
                              "region4" = "Countryside",
                              "cntryIT" = "Italy",
                              "cntryPL" = "Poland"),
       groups = list("Crisis specific" = 1:4,
                     "Socio-demographics" = 5:9,
                     "Income (ref: comfortable)" = 10:12,
                     "Region (ref: big city)" = 13:15,
                     "Country FE (ref: Germany)" = 16:17),
       include.deviance = FALSE,
       include.loglik = FALSE,
       # file=paste0(fig_path, "general_participation.tex"),
      custom.note	= "%stars")

```

Results are based on the first wave, weighted by the socio-demographic weight. Let's see the figure with the odds ratios:

\newpage

### Figure 2, Appendix A

```{r, include = F, echo = F, message = F, warning=FALSE}

mod1_res <- tidy(M1)
mod1_res <- mod1_res %>% 
  mutate(type="Not contacted")
mod1_ci <- as.data.frame(confint(M1))
mod1_ci$term <- rownames(mod1_ci)
mod1_res <- merge(mod1_ci, mod1_res, by="term", all=TRUE)
mod1_ci <- as.data.frame(confint(M1, level=0.9))
mod1_ci$term <- rownames(mod1_ci)
mod1_res <- merge(mod1_ci, mod1_res, by="term", all=TRUE)

mod2_res <- tidy(M2)
mod2_res <- mod2_res %>% 
  mutate(type="Contacted")
mod2_ci <- as.data.frame(confint(M2))
mod2_ci$term <- rownames(mod2_ci)
mod2_res <- merge(mod2_ci, mod2_res, by="term", all=TRUE)
mod2_ci <- as.data.frame(confint(M2, level=0.9))
mod2_ci$term <- rownames(mod2_ci)
mod2_res <- merge(mod2_ci, mod2_res, by="term", all=TRUE)

probs <- bind_rows(mod1_res, mod2_res)

```


```{r, include = T, echo = F, message = F, warning=FALSE, fig.width=5, fig.height=5}

plot_dat <- probs %>% 
  filter(term!="(Intercept)") %>% 
  filter(!grepl("cntry", term)) %>% 
  mutate(odds = exp(estimate),
         odds.outter.upper = exp(`2.5 %`),
         odds.outter.lower = exp(`97.5 %`),
         odds.inner.upper = exp(`5 %`),
         odds.inner.lower = exp(`95 %`)) %>% 
  mutate(term=case_when(term=="edu_lev21" ~ "High educ.",
                        term=="gndr2" ~ "Female",
                        term=="region2" ~ "Mid. town",
                        term=="region3" ~ "Small town",
                        term=="region4" ~ "Countryside",
                        term=="org_membership1" ~ "Memb. in civ. soc. org",
                        term=="age" ~ "Age",
                        term=="age2" ~ "Sq(Age)",
                        term=="infection1" ~ "Infection (self/ in household)",
                        term=="home_child1" ~ "Children at home",
                        term=="home_care1" ~ "Care work at home",
                        term=="neg_income_chng" ~ "Negative income change",
                        term=="income2" ~ "Coping",
                        term=="income3" ~ "Difficult",
                        term=="income4" ~ "Very difficult",
                        TRUE ~ term)) 

titles <- plot_dat %>%
  filter(term %in% c("Age", "Coping", "Difficult", "Female")) %>%  # I just need some unique values
  mutate_at(vars(c(2:9), c(11:15)), ~ NA)

titles$term <- rep(c("Crisis specific", "Socio-demographics", "Income (ref: comfortable)", "Region (ref: big city)"), 1)

plot_dat <- bind_rows(plot_dat, titles) %>% 
  mutate(term=factor(term, levels = c("Crisis specific",
                                      "Infection (self/ in household)",
                                      "Children at home",
                                      "Care work at home",
                                      "Negative income change",
                                      "Socio-demographics",
                                      "Memb. in civ. soc. org",
                                      "Age", "Sq(Age)", "High educ.", 
                                      "Female", 
                                      "Income (ref: comfortable)",
                                      "Coping", "Difficult", "Very difficult",
                                      "Region (ref: big city)",
                                      "Mid. town", "Small town", "Countryside"))) %>%
  mutate(term=fct_rev(term)) %>% 
  mutate(new_labels = as.character(term)) %>% 
  mutate(DV_name = "Volunteering") 

space_between_bars <- 0.5

res1 <- ggplot(plot_dat, aes(y=term, color=type)) +
  geom_vline(xintercept = 1, linetype="dashed", color="gray60") +
  geom_point(aes(x=odds, shape=type), position=position_dodge(width=space_between_bars), size=2) +
  geom_linerange(aes(xmin=odds.outter.lower, xmax=odds.outter.upper, 
                     group=type), size=0.5, position=position_dodge(width=space_between_bars), key_glyph = "path") +
  geom_linerange(aes(xmin=odds.inner.lower, xmax=odds.inner.upper, 
                     group=type), size=1.1, position=position_dodge(width=space_between_bars), key_glyph = "path") +
  xlab("Odds Ratio") +
  scale_color_manual(name="", values=c("#758fc7", "#cc592d"), guide=guide_legend(reverse = TRUE, nrow=1, byrow=TRUE)) +
  scale_shape_discrete(name="", guide=guide_legend(reverse = TRUE, nrow=1, byrow=TRUE)) +
  scale_y_discrete(labels = function(y) stringr::str_wrap(y, width = 22)) +
  theme_linedraw() +
  theme(legend.title=element_blank(), 
        legend.position="bottom",
        legend.key.width = unit(1,"cm"),
        axis.title.y = element_blank(),
        axis.text.y=element_text(face=c(rep('plain', 3), 'bold', 
                                        rep('plain', 3), 'bold', 
                                        rep('plain', 5), 'bold', 
                                        rep('plain', 4), 'bold'),
                                 hjust=1),
        # axis.text.y = element_text(face = c(rep(2, 4))),
        axis.ticks.y = element_blank())

# ggsave(plot=last_plot(),
#        filename="volunteering_by_contact.png",
#        path = fig_path,
#        width=5,
#        height=5,
#        scale=3,
#        dpi=400,
#        units="cm")

res1 <- res1 +
  facet_wrap(~DV_name)

res1

```

